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INTRODUCTION 


I . 


This  is  the  seventh  semi-annual  report  and  covers  work 
performed  from  October  1964  through  March  1965.  During  this 
period,  an  amendment  was  negotiated  to  extend  contract  AF33(657)- 
12447  from  17  February  1965  through  17  February  1966.  Under 
terms  of  this  amendment,  semi-annual  reports  will  be  replaced 
by  quarterly  reports;  the  present  report  combines  the  semi¬ 
annual  previously  required  and  the  first  quarterly  report 
required  under  the  amended  contract. 

Specific  analyses  made,  for  which  results  have  been 
reported  to  AFTAC ,  are  discussed  in  Section  II  under  descriptive 
headings;  work  currently  in  progress,  which  has  not  reached 
the  stage  where  results  are  available,  is  discussed  in  Section 
IV,  Work  previously  completed  is  mentioned  only  as  it  relates 
to  analyses  in  progress  during  this  reporting  period. 

II .  WORK  COMPLETED 

A .  Analysis  of  Signal  Variance 

The  purpose  of  the  present  work  has  been  to  give  a 
technical  evaluation  of  the  analysis  of  variance  as  a  method  of 
automatic  signal  detection.  The  analysis  of  variance  has  been 
programmed  for  the  digital  computer  and  the  resulting  program 
may  be  considered  an  array  processor  for  the  purpose  f  signal 
detection  only. 

The  analysis  of  variance  considered  as  a  method  of 
signal  detection  by  Melton  and  Bailey^  depends  upon  the  simi¬ 
larity  of  signals  across  channels,  as  compared  to  the  noise, 
for  array  recorded  data.  For  weak  signals,  prefiltering  is 
necessary  as  the  similarity  of  the  long-period  noise  across 
the  array  dominates  the  results.  Thus  , the  computational  task 
is  to  filter  out  the  long-period  noise,  hopefully  leaving  only 
incoherent  noise  across  channels.  Then  the  analysis  of  variance 
computation  provides  an  output  function  of  time  which  measures 
the  likelihood  of  a  signal  being  present  in  the  noise. 

One  important  feature  of  the  above  detection  technique 
is  a  lack  of  assumptions  about  the  signal.  Most  of  the  classical 
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Melton,  Een  S.  and  Bailey,  Leslie  F.,  'Multiple  Signal 
Correlators,"  Geophysics ,  Vol .  XXII,  No.  3,  July  1957. 


work^  in  optimum  filtering  or  detection  theory  assumes  that 
both  signal  and  noise  are  stochastic  processes  with  known 
spectra.  The  resulting  optimum  processing  was  given  by 
Weiner.  If  one  assumes  the  signal  is  a  fixed  known  wavelet, 
then  the  optimum  filter  is  called  a  matched  filter.  Modi¬ 
fications  of  the  above  techniques  have  been  studies,  such 
as  assuming  that  the  signal  is  repeating  in  the  noise  or  is 
known  apart  from  a  set  of  parameters  which  have  to  be  estimated 
from  the  data.  None  of  these  techniques  provides  a  suitable 
basic  assumption  regarding  seismic  signals.  If  the  signal 
trace  is  regarded  as  a  stochastic  process  then  it  is  non- 
itationary  and  zero  with  high  probability.  If  the  signal 
trace  is  considered  deterministic,  then  the  variation  due  to 
recording  station,  epicenter  location,  or  type  of  event  is 
difficult  to  model  with  a  relatively  small  number  of  para¬ 
meters.  Consequently,  it  is  considered  desirable  to  examine 
..hose  techniques  which  make  no  assumptions  at  all  regarding 
the  signal. 


The  logical  steps  of  computation  for  the  method  of 
detection  considered  in  this  study  were  subsetting  data, 
signal  addition,  filtering,  and  analysis  of  variance  computa¬ 
tion.  These  steps  have  been  programmed  individually  and 
combined  into  a  single  program  by  a  system  called  "FLAP". 


The  basic  operation  the  program  must  repeatedly 
make  is  to  compute  Fisher's  F  statistic  for  testing  the 
hypothesis  that  no  column  effects  (no  signal)  exist  in  a 
2-way  table  with  one  observation  per  cell  is  given  by: 


,s(r-l) 
P  r(s-l) 
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1  3 


Y< 
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ij 
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where  y.  .  is  the  observed  value  for  the  i  *  channel,  and  the 
jth  sampling  time,  r  the  number  of  channels,  and  s  the  number 
of  time  increments  at  which  samples  are  measured. 
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Pergamon  Press,  1960. 


Several  parameters  such  as  sampling  rate,  bard 
pass  filtering,  and  signal-to-noise  ratio  were  tested  for 
effect  on  the  computation  of  F.  From  the  following  table, 
the  effect  of  changing  the  sampling  rate  is  seen  to  be 
minimal . 


DATA  SAMPLE  (Filtered  TFO  Noise) 


SAMPLING  RATE 
(Samples  Per  Second) 

1 

2 

3 

4 

100 

1.94 

2.16 

2  .24 

1.35 

50 

1.94 

2.17 

2.24 

1.35 

20 

1.95 

2.19 

2  .25 

1.35 

10 

1.98 

2  .21 

2.27 

1.35 

DATA  SAMPLE  (Filtered  TFO  Signal  and  Noise) 


1 

2 

3 

4 

100 

47.31 

67.58 

36.60 

43  .l'i 

50 

47  .40 

67.70 

36.70 

43.30 

20 

47.65 

68.13 

36.87 

43.49 

10 

48.51 

68.88 

37.23 

43.78 

Effect  of  Sampling  Rate  on  F 


It  is  clear  that  100  samples  per  second  is  unnecessary.  In 
fact,  the  results  indicate  that  10  samples  per  second  is 
sufficiently  accurate  and  the  effective  sampling  rate  was 
frozen  at  this  value  for  the  remainder  of  the  results. 

From  results  of  applying  the  detection  procedure 
to  raw  data  and  band-pass  filtered  data,  it  was  obvious  that 
filtering  is  a  critical  part  of  the  method.  The  band-pass 
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filter  with  cut-offs  at  0.75  cps  and  4  cps  was  selected  as 
standard  for  the  present  report.  Similarly,  several  values 
of  window  width  were  tried  with  the  expected  result  that 
the  optimum  window  width  is  the  length  of  the  signal. 

Figure  1  illustrates  the  effect  of  S/N  ratios  on 
F  computations  with  signals  of  varying  complexity  at  several 
stations . 


The  F's  were  calculated  using  5  second  windows 
which  included  the  array  stations  first  second  of  signal 
plus  4  seconds  of  preceding  noise.  In  general  the  slope 
of  the  curves  for  simpler  signals  is  greater,  but,  in  any 
case,  the  expected  contrast  in  F  values  can  be  expected  to 
be  appreciably  greater  than  the  S/N  ratio  itself.  The 
background  noise  given  cn  the  figure  is  an  average  of  23 
values  of  F  for  the  noise  preceding  the  signal. 

Eighteen  Leleseismic  signals  were  selected  on  the 
basis  of  small  signal-to-noise  ratios.  No  velocity  filtering 
was  done  before  processing.  Three  examples  of  the  results 
are  shown  in  Figure  2.  In  each  case,  the  dotted  trace  gives 
the  F  values  and  the  solid  trace  the  unphased  sum  of  all 
seismograms  from  the  array.  Timing  lines  are  spaced  at 
5-second  intervals . 

In  an  automatic  detection  system,  one  would  have 
to  assume  sore  critical  value  F  above  which  the  F-trace 
would  sound  an  alarm.  In  this  experiment,  the  theoretical 
false  alarm  rate  was  surpassed  somewhat.  This  was  to  be 
expected  as  we  knew  the  band-pass  filtering  would  not  com¬ 
pletely  eliminate  the  correlations  in  time,  and  yet  our 
assumptions  were  that  the  observations  in  time  are  uncorre¬ 
lated.  One  should  be  able  to  adjust  the  critical  values  to 
take  into  consideration  this  correlation.  This  problem  has 
been  consiaered  in  the  literature  by  G.  E.  P.  Box3  for 
correlations  extended  only  one  unit  and  some  calculations 
were  necessary  to  extend  this  theory  to  the  present  problem. 

We  also  studied  the  effects  of  velocity  filtering 
the  data  before  computing  F  values. 
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SIGNAL  TO  NOISE  RATIO 


Increase  In  F  As  A  Function  Of  S/N 


Figure  1 


Stumog'om  No.  2728 
UBO 


F  Values  and  Corresponding  Array  Sums  for  Teleseismic  Signals 


The  results  indicated  that  time-shifting  the  data 
is  not  essential  for  the  detection  procedure,  although  signifi¬ 
cant  improvement  may  be  possible  for  some  signals  at  low 
apparent  velocities.  The  cost  of  this  improvement  is  that  F 
values  for  a  grid  of  velocities  and  epicentral  azimuth  would 
have  to  be  computed  for  each  window  of  data.  More  work  is 
needed  to  ascertain  whether  the  gain  will  justify  this  cost, 
although  in  some  cases  when  an  operational  array  is  to  be 
beamed  at  a  single  epicentral  region,  the  additional  cost 
would  be  minimized. 

We  have  concluded  fhat  the  method  of  signal  detection 
examined  is  very  effective  when  the  background  noise  is  incoher¬ 
ent  in  the  frequency  band  of  interest .  The  background  noise 
seemed  to  satisfy  this  criterion  with  the  exception  of  TFO 
where  the  average  background  values  of  F  were  at  least  two 
times  the  average  values  at  the  other  arrays.  This  may  be 
due  to  a  completely  different  character  of  the  noise  field 
at  TFO. 


There  is  evidence  that  very  accurate  shifting  would 
improve  the  significance  of  the  F  values.  Not  only  are  the  F 
values  increased  due  to  shifting  when  signals  are  present, 
but  the  background  noise  values  are  systematically  reduced  by 
shifting  in  the  absence  of  signal.  The  present  program  com¬ 
putes  the  shift  in  units  of  0.05  seconds  and  truncates  to  the 
nearest  integer.  Future  programs  should  be  written  to  inter¬ 
polate  as  an  option  to  allow  a  more  accurate  view  of  the  effect 
of  time  shifting. 

The  large  increase  in  S/N  for  the  F  "seismogram"  has 
been  seen  by  comparing  the  peak  values  of  F  to  the  mean  back¬ 
ground  values.  The  parameter  x^  was  also  studied  for  com¬ 
parison.  Although  the  detectors  F  and  were  both  effective 
as  detectors  for  the  data  examined,  the  F  peak  values  ex¬ 
ceeded  the  mean  background  values  by  a  consistantly  larger 
number  of  standard  deviations  than  did  the  X^ •  The  F  criteria 
is  amplitude -dependent  where  the  X*^  values  are  not.  Thus,  one 
would  expect  the  F  detector  to  be  better  for  strong  signals. 

The  data  show  the  F  values  are  better  for  weak  signals  as  well. 
Unless  there  are  other  compelling  reasons,  such  as  cost  of 
computation  or  possiblity  of  on-line  computation  of  X^  and  not 
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F,  the  results  of  this  report  confirm  the  theoretical  result 
that  F  should  be  an  optimum  detection  criterion. 

The  following  work  is  recommended  if  the  F  method  is 
to  be  considered  further  for  noise  elimination  with  the  LASA 
data: 


1.  Re-program  the  analysis  of  variance  making  the 
following  changes: 

A.  The  filtering  should  be  made  an  integral 
part  of  the  analysis  of  variance  computation. 

B.  The  program  should  handle  a  maximum  of  525 
channels  with  the  computation  being  done 
"in  time",  i.e.,  the  calculations  necessary 
for  a  vector  of  observations  at  a  given  point 
in  time  should  be  completed  at  once,  then 
the  next  vector  in  time  processed,  and  so  on. 

C.  Shifting  with  interpolation  should  be 
included  as  an  option. 

D.  The  program  should  be  written  for  use  on  the 
1604  but  with  the  idea  in  mind  that  the  same 
flow  charts  and  analysis  would  be  suitable, 
with  minor  changes,  for  programming  the  on¬ 
line  system  at  LASA. 

2.  The  program  described  in  (1)  should  be  used  to 
investigate  sensitivity  to  velocity  phasing  on 
larger  arrays. 

3.  The  existing  program  should  be  used  to  study 
various  filters  for  the  different  arrays.  Optimum 
single  channel  filtering  in  frequency  could  be 
used  and  the  results  compared  with  those  of  the 
present  report.  (Optimum  multi-channel  filtering 
may  introduce  correlations  across  channels.) 

B .  Solutions  of  the  General  Wave  Equation 

Midway  in  the  current  reporting  pen-  formal  work  on 
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this  project  was  suspended  indefinitely  in  favor  of  studies 
that  bear  promise  of  more  immediate  benefit  to  the  detection- 
identification  problem.  However,  several  significant  results 
were  reported  before  the  work  was  terminated. 

Green's  function  for  the  reduced  wave  equation 
(Helmholtz  equation)  in  a  spherical  annular  domain  with  Dirichlet 
boundary  conditions  was  derived.  The  convergence  of  the  series 
solution,  representing  Green's  function,  was  then  established. 
Finally,  it  was  shown  that  Green's  function  of  the  Dirichlet 
problem  reduces  to  Green's  function  for  the  exterior  of  a  sphere 
as  given  by  Franz  and  Etienne^,  when  the  outer  radius  is  moved 
towards  infinity,  and  when  a  special  position  of  the  cooredinate 
system  is  chosen . 

Green's  function  has  been  constructed  here  in 
reference  to  an  asymmetrically  located  coordinate  system, 
different  from  Franz's  and  Etienne's  results  where  an  axial 
symmetric  position  was  employed.  The  reason  for  this  general¬ 
ization  lies  in  the  fact  that  in  contrast  to  axial  symmetric 
systems,  disturbance  fields  due  to  earthquakes,  for  example, 
are  most  probably  not  axially  symmetric  relative  to  the  approxi¬ 
mate  spherical  shape  of  the  earth.  This  type  of  Green's  function 
permits  construction  of  any  general  asymmetric  discrete  source 
system  inside  a  spherical  layer  Therefore,  a  comparison  of  at 
least  pure  steady  state  body  waves  and  their  scattering  cha:  ac- 
teristics  inside  fluid  sphere  layers  with  free  boundaries, 
where  the  waves  are  created  from  point  symmetric  or  axially 
symmetric  sources  on  the  one  hand,  and  asymmetric  sources  on  the 
other  hand,  can  now  be  made.  According  to  an  idea  advanced  by 
Jeffreys,  this  comparison  indicates  the  body  wave  scattering 
characteristics  occunng  in  elastic  spheres. 

The  corresponding  solutions  and  proofs  were  given 
for  a  circular  annular  domain  with  the  Dirichlet,  the  radiation 
(mixed),  and  Neumann  boundary  conditions  being  derived.  The 
convergence  of  the  series  representing  Green's  functions  was 
then  established.  Finally  it  was  shown  that  these  functions 
reduce  to  Green's  function  for  the  exterior  of  a  circle  as 
g.  'en  by  Franz  and  Etienne,  when  the  outer  radius  is  moved 
towards  infinity. 
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Currently,  in  seismology,  numerical  analysis  is 
applied  to  obtain  normal  mode  solutions  for  propagating  waves 
in  a  plane,  parallel  homogeneously  layered  structure.  The 
solutions  are  obtained  rapidly  and  with  high  accuracy.  More 
recently,  normal  mode  solutions  in  a  spherical  earth  model 
have  been  obtained  numerically  for  spheres  composed  of  con¬ 
centric  shells.  However,  in  this  case,  there  was  considerable 
difficulty  in  obtaining  solutions  rapidly  and  accurately.  The 
main  difficulties  were  posed  by  the  requirement  of  an  excessive 
number  of  shells  to  accurately  represent  the  whole  interior  of 
the  earth,  and  computational  difficulties  in  evaluating  Bessel 
functions  with  sufficient  accuracy,  especially  at  high  frequen¬ 
cies.  We  conducted  an  analyis  to  explore  two  avenues  toward 
alleviating  these  difficulties.  In  particular  we: 

1.  Represented  the  elastic  parameters  in  each  shell 
as  continuously  variable  quantities  in  order  to 
greatly  reduce  the  number  of  shells  required  to 
accurately  represent  the  interior  of  the  earth. 

2.  Rigorously  evaluated  a  transformation  through  which 

a  solution  in  plane  parallel  layers  could  be  utilized 
to  approximate  the  solution  for  spherical  shells. 

In  the  latter  case,  a  special  kind  of  inhomogeneous 
medium  has  been  analyzed  in  which  the  reciprocal  radii  trans¬ 
formation  for  spherical  inversion  is  invariant  with  respect 
to  the  reduced  wave  equation  describing  the  motion  in  a  sphere. 
Thus,  by  transforming  a  shell  from  a  region  inside  of  a  large 
inversion  sphere  to  a  region  outside  that  sphere,  the  boundary 
conditions  can  be  accurately  represented  by  those  for  plane 

if  pseudo  amplitudes  are 
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parallel  layers.  In  particular, 
defined  as 


-z 


C'l  =  C1  (  2a 


•) 


-z 


C '  2  =  c2  exp  (  2a 


and  a  pseudo  velocity  is  defined 


P0r|  ”  2  2 
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where  a  is  the  radius  of  the  inversion  sphere,  v  either  the  shear 
or  compressional  propagation  velocity,  p  the  density,  and  uu  the 
angular  frequency,  the  wave  potential  is  reduced  exactly  to  the 
form  of  a  plane  wave  potential  in  a  flat  layered  structure,  as 
follows : 


R( z ' )  =  C'  exp 


-1  z' J  ,  (0  <  z'  <  ®) 

Haskell's5  exact  method  for  plane  parallel  layers 
can  now  be  applied  in  principle  to  obtain  the  approximate 
solution  of  a  layered  elastic  sphere.  To  obtain  free  periods 
of  oscillation  of  the  layered  sphere  we  take  a^k2  =  n(n  +  l), 
which  for  each  wave  mode  specified  by  Ci  (id)  leads  to  eigen- 

Ci(uj)2  2 

frequencies  given  by  roots  of  the  equation  n(n+l)  ~  -  w  -  0. 

a 

The  completeness  of  the  exact  solutions  obtained  using  Haskell's 
method  will  depend  upon  the  treatment  of  the  boundary  condition 
at  infinity,  which  corresponds  to  r  =  0  in  the  region  0  <  r  <  a. 

It  is  important  to  remember  that  Sommer f eld  showed 
that  spherical  inversion  cannot  be  applied  to  the  usual  Helmholtz 
equation.  Here  it  has  been  shown  that  spherical  inversion  can  be 
applied  to  at  least  one  special  case  of  the  reduced  wave  equation 
with  variable  elastic  parameters. 

In  order  to  determine  the  effect  of  diffraction 
around  the  earth's  core,  we  studied  SH  waves  in  a  model  con¬ 
sisting  of  a  single  spherical  elastic  layer  which  surrounds  a 
liquid  core.  We  considered  two  SH  sources  located  inside  the 
elastic  layer,  such  that  the  field,  due  to  the  sources,  is 
axial  ‘'-vmmetric  and  steady  state.  Thus,  the  sources  have  been 
assu  in  a  form  which  assures  equilibrium  of  the  elastic 
fluid  system  at  all  times.  The  effect  of  gravity  has  been 
neglected . 
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We  commenced  by  deriving  the  exact  solution  due  to  a 
single  SH  source  in  form  of  Green's  function  and  investigating 
the  denominator  D  of  this  series  solution  in  relation  to  the 
complex  Watson  integral.  The  consideration  of  a  single  SH  source 
was  here  sxifficient  as  all  other  SH  source  distributions  were 
immediately  obtained  by  superposition.  For  this  purpose,  it  was 
necessary  to  consider  the  Green's  function  as  a  function  of  the 
index,  which  now  appears  as  an  independent  variable,  and  to  study 
those  indices  of  the  functions  appearing  ir.  D  for  which  D  vanishes. 
Calling  the  index  of  the  functions  appearing  in  D  "v",  it  was 
demonstrated  that  the  v- zeros  of  D  were  similar  to  the  ones  obtained 
for  the  Green's  function  for  the  reduced  wave  equation  in  an 
annular  domain  and  in  a  spherical  layer  with  Dirichlet's,  Neumann's, 
or  mixed  boundary  conditions.  The  v-zeros  were  found  to  be  either 
real  or  complex.  Futhermore,  for  real  arguments  of  the  functions 
appearing  in  D ,  it  was  established that  D  had  a  countable  number 
of  v-zeros  of  which  only  at  most  a  finite  number  were  real.  All 
of  these  v-zeros  were  shown  to  be  simple,  except  for  that  which 
occurred  at  =  -  h,  which  can  appear  only  as  a  zero  of  the 

second  order.  Asymptotic  expressions  were  then  derived  for  the 
two  dimensional  counterpart  for  large  imaginary  v-zeros. 


In  addition,  it  was  proven  that  there  are  no  real  v- 
zeros  for  arbitrarily  large  modulus.  Also,  D  was  found  to  be 
neither  an  odd  nor  even  function  of  v  for  the  general  case;  however, 
symmetry  of  the  v-zeros  with  respect  to  the  real  v-axis,  of  course, 
does  exist.  The  conclusion  drawn  from  these  studies  on  the  v-zeros 
indicates  that  the  Watson  transformation  is  only  directly  applicable 
in  those  physical  tuations  where  there  are  no  real  v-zeros.  For 
these  special  cases  a  core  comparison  is  then  carried  out,  whereby 
the  two-dimensional  counterpart  to  the  real  three-dimensional 
elasticity  problem  has  been  chosen  for  reasons  of  simplicity. 

These  studies  reveal  that  the  only  difference  for  certain  simple 
core  changes,  namely  =  0  compared  to  0,  is  essentially 

a  shift  in  the  pure  imaginary  v-zeros  along  that  axis.  In  this 
instance,  the  distribution  of  the  imaginary  v_Zeros  is  hardly 
effected . 


In  comparison,  in  problems  in  seismology,  for  the  type 
of  layer  and  core,  and  nature  of  source  considered  here,  the  waves 
are  unattenuated.  Also,  due  to  the  periods,  wave  numbers,  etc.,  of 
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rhe  waves  which  occur  in  the  earth,  and  the  dimensions  of  the 
layers,  only  real  v -zeros  exist  here.  In  this  situation,  like 
in  the  one  considered  above,  where  there  are  only  complex  (in 
three  dimensions)  or  purely  imaginary  (in  two  dimensions)  v -zeros, 
the  elastic  layer  acts  as  a  waveguide.  This  is  due  to  the  fact 
that  there  are  •  o  shear  waves  transmitted  to  che  outside  cf  the 
layer  nor  into  the  core.  In  these  cases,  the  Watson  transformation 
is  no  longer  direct lv  applicable  to  the  complete  solution.  A 
decomposition  of  the  complete  solution  is  therefore  carried  out 
on  similar  lines  as  was  done  by  Van  der  Pol  and  Bremmer^  by 
constructing  the  solution  anew  in  terms  of  infinite  sums  of 
reflections  inside  the  layer  with  the  aid  of  reflection  coef¬ 
ficients.  These  coefficients  are  consequently  computed  and  it 
is  shown  that  Watson's  transformation  can  now  be  again  applied 
to  each  term  of  the  new  series,  as  the  individual  terms  of  the 
new  series  are  finally  free  of  real  v-zeros.  In  order  to  facili¬ 
tate  the  application  of  the  Watson  transformation  to  expressions 
with  multiple  v-2eros,  which  characterize  the  terms  of  the  new 
series,  the  necessary  formula  based  on  Cauchy's  theorem  and  the 
Laurent  expansion  is  derived. 

In  conclusion,  a  concise  analysis  of  the  Cauchy 
formula  and  the  Laurent  expansion  yielding  the  residues  has 
been  given,  as  the  residues  in  a  layer  problem  are  much  more 
complicated  than  the  residues  in  problems  with  a  single  surface 
boundary.  The  method  derived  in  this  work  brings  the  problem 
to  a  stage  where  numerical  computations  can  now  be  contemplated. 
This  is  possible  because  each  of  the  denominators  of  the  geo¬ 
metric  series  terms  has  now  v-zeros  whose  analysis  has  been 
performed  earlier  in  other  investigations  (Duvalo  -  Jacobs)^. 

With  this  result  it  is  finally  possible  to  separate  approxi¬ 
mately  surface  (diffraction)  waves  from  reflected  waves  and 
determine  their  nature  and  characteristic  features. 

C .  Model  Studies 

We  have  an  operational  program  to  model  the  effect  of 
reverberation  of  teleseismic  signals  in  the  crust,  both  at  receiver 
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and  the  source.  Although  we  are  no  longer  doing  development 
work,  we  have  used  this  program  to  demonstrate  some  simple, 
albeit  general,  conclusions  concerning  teleseismic  signals. 

The  wave  forms  in  Figure  3  were  plotted  to  furnish  an 
elementary  insight  into  the  effect  of  the  free  surface  on  the 
primary  wave  from  a  shallow  explosion. 

The  wave  forms  shewn  in  the  upper  half  of  Figure  3 
are  simply  the  sum  of  two  identical  wavelets,  one  of  which  is 
inverted  and  delayed  by  the  intervals  shown  (At).  The  results 
are  interpreted  to  be  simulated  downward  travelling  F -pulses, 
perturbed  by  the  surface  reflected  pP-pulses.  Inasmuch  as  the 
apparent  period  of  the  basic  wavelet  is  one  second,  the  delay 
times  may  be  interpreted  as  the  two-way  travel  paths  from  the 
explosion  to  the  earth's  surface,  expressed  as  fractions  of 
wavelengths . 

An  analogous  set  of  signals  is  shown  in  the  lower 
half  of  Figure  3,  in  which  the  two  signals  are  added  in  phas* 
and  delayed  as  before.  These  wave  forms  may  be  interpreted 
as  incoming  teleseisms,  originally  of  the  same  simple  pulse 
shape,  recorded  at  depth  below  the  free  surface. 

Although  the  wave  shape  varies  considerably  as  a 
function  of  depth,  there  is  no  quantity  or  characteristic  which 
can  be  measured  and  related  to  depth.  The  only  way  to  relate 
wave  shape  to  depth  for  the  comparatively  shallow  events 
depicted  in  this  exercise  would  be  to  match  wave  forms.  This 
technique  also  would  probably  be  unsuccessful  in  practice 
due  to  complications  beyond  our  control. 

In  a  similar  exercise,  we  have  studied  the  possible 
effects  of  reverberation  on  signal  amplitude.  When  the  signal 
amplitude  recorded  at  each  of  several  stations  for  a  given 
event  is  plotted  as  a  function  of  the  epicentral  distance,  there 
is  considerable  scatter  in  the  data.  The  same  is  true  even  when 
special  care  is  taken  to  insure  that  the  points  represent  the 
arrivals  of  the  same  single  phase.  It  was  the  purpose  of  this 
study  to  investigate  in  a  very  simple  way  what  possible  variation 
in  signal  amplitude  we  might  reasonably  expect  from  changes  in 
the  near  surface  rocks  and,  in  particular,  from  changes  in  the 
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sedimentary  section  from  station  to  station 


The  basic  contributing  unit  in  the  geologic  section 
is  obviously  the  single  stratigraphic  layer  and  we  are  intuitively 
led  to  believe  that  the  most  important  layer  is  the  surface 
layer  where  the  topside  reflection  coefficient  is  unity  or,  in 
practice  at  least,  near  unity.  Therefore,  we  commenced  our 
study  with  a  single  uniform  layer.  The  density  within  the 
layer  is  2.35  grams  per  cubic  centimeter  and  the  velocity  of 
compress ional  waves  is  3.1  km  per  second.  In  the  underlying 
material  the  corresponding  density  and  velocity  are  2.60  and 
6.03 . 


Figure  4  shows  the  effect  on  the  amplitude  of  the 
signal  at  the  surface;  the  amplitudes  are  normalized  so  that  the 
amplitude  for  zero  thichness  is  unity  and  the  thicknesses  are 
measured  in  terms  of  wavelength  X.  We  could  have  predicted 
that  we  would  have  obtained  a  maximum  for  a  resonant  thickness 
equal  to  quarter-wavelengta  and  this  curve  indicates  how  large 
the  maximum  amplitude  is.  For  velocities  shown  in  Model  1, 
and  for  teleseisms  peaking  at  one  cycle  per  second,  the  quarter- 
wr-.velength  resonant  thickness  is  775  meters  and  the  signal  amp¬ 
litude  is  up  91%  over  its  level,  were  tl  top  layer  completely 
removed,  it  should  be  noted  that  tho  amplitude  in  question  is 
the  maximum  peak-to-peak  amplitude  which  is  usually  measured 
from  the  first  peak  to  the  second  trough. 

Having  verified  that  the  maximum  amplitude  occurs  at 
a  quarter-wavelength  thickness,  it  is  now  useful  to  examine  how 
the  amplitude  varies  when  the  acoustic  impedance  is  changed. 

To  gain  insight  into  this  problem,  we  repeated  the  first  exer¬ 
cise,  simply  changing  the  density  and  compressional  velocity  in 
the  upper  layer  to  2.15  and  2.0.  These  are  about  the  lower 
limits  of  the  parameters  that  are  found  for  a  reasonably  thick 
surface  layer  in  nature.  The  results  are  also  shown  in  Figure 
4.  As  in  the  first  case,  the  maximum  amplitude  occurs  at  a 
thickness  equal  to  the  quarter-wavelength.  However,  the  amp¬ 
litude  at  the  maximum  has  changed,  as  has  the  amplitude  for 
thicknesses  of  the  order  of  a  wavelength  and  larger. 

The  limiting  amplitude  for  large  thicknesses  is 
obviously  governed  by  the  acoustic  impedance;  the  ratio  of  the 
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Signal  Amplitude  As  A  Function  Of  Thickness  Of  A  Single  Surface  Layer 


two  limiting  amplitudes  is  1.30  which  is  equal  to  the  square 
root  of  the  inverse  ratio  of  the  corresponding  acoustic  impedance 
The  amplitude  at  resonance  depends  not  only  on  the  acoustic  imped 
ance  but  als''  on  the  reflection  coefficient  at  the  boundary. 

Thus,  the  resonant  amplitude  for  Model  2  is  1.62  times  the 
resonant  amplitude  for  Model  1. 

We  experimented  with  several  multi-layered  models, 
including  those  with  velocity  inversions,  but  could  find  none 
that  would  give  us  appreciably  larger  amplitude  effects.  In 
conclusion,  we  have  found  that  resonance  is  the  most  important 
factor  in  signal  amplitude  variation  resulting  from  changes  in 
near  surface  geology.  To  a  lesser  extent,  the  acoustic  imped¬ 
ance  of  the  surface  rock  is  also  important,  a  smaller  acoustic 
impedance  at  the  surface  leading  to  a  large!  signal  amplitude. 

As  an  order  of  magnitude,  the  signal  amplitude  may  vary  by  a 
factor  of  three  or  four  due  to  such  causes  in  an  idealized 
case;  in  practice,  the  variation  is  probably  less. 

In  support  of  current  interest  in  the  performance  at 
various  sites,  we  used  the  above  technique  in  predicting  the 
comparative  amplitudes  of  teleseismic  signals  received  at  a 
number  of  stations.  These  computations  varied  from  the  previous 
ones,  however,  in  that  the  reference  signal  of  unit  amplitude 
was  now  incident  on  the  Moho  and  not  the  top  of  the  granite 
layer.  In  other  respects,  the  present  and  previous  computations 
were  the  same. 


The  models  and  results  are  shown  in  Figure  5.  The 
models  are  based  on  as  precise  data  as  is  available,  plus,  in 
some  cases,  some  guidance  from  a  knowledge  of  regional  geology. 
In  each  of  the  models  shown,  the  surface  layer  is  sedimentary 
with  a  distinctly  low  velocity.  In  areas  of  igneous  outcrops, 
such  as  those  classically  chosen  as  ideal  seismometer  sites,  the 
models  would  be  essentially  as  those  shown  without  the  sedi¬ 
mentary  cover  of  velocities  below  6  km  per  second.  The  thick¬ 
ness  of  the  surface  layer  would  be  such  that  there  would  be 
insignificant  resonance. 

In  each  case,  the  amplitude  of  the  signal  incident 
on  the  Moho  from  below  are  considered  to  be  unity.  This  means 
that,  if  there  were  no  overlying  crust,  the  amplitude  would  be 
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Figure  5 


2.0  at  the  free  surface.  The  corresponding  amplitudes  at  the 
four  sites  are  3.64  for  Las  Cruc  ,s,  N.  M.  (LC  NM>)  ,  2.40  for 
Delhi,  N.  Y.  (DH  NY) ,  3.88  for  Winner,  C .  D.  (WN  SD)  ,  and  3.00 
for  Uinta  Basin  Observatory,  Utah  (UBO) . 

If  there  were  no  sedimentary  layer  and  the  surface 
layer  were  granite  with  a  velocity  of  6.03  and  density  2.60,  we 
used  the  acoustic  impedances  to  get  2.64  as  the  corresponding 
amplitude.  In  computing  this  figure,  we  assumed  no  loss  due 
to  reflections  at  intermediate  surfaces;  thus  the  amplitude  of 
2.64  must  be  considered  a  reasonable  approximation  but  still  an 
upper  limit.  The  thicknesses  in  this  crust  are  such  that  we 
can  disregard  reverberations  in  considering  *-he  amplitude. 

If  we  now  wish  to  normalize  the  amplitudes  in  the 
sedimentary  models  to  that  in  the  shield  model  taken  as  unity, 
we  have  1.38  for  LC  NM,  0.91  for  DH  NY,  1.47  for  WN  SD,  and 
1.14  for  UBO. 

In  general,  we  can  conclude  that  the  signals  ampli¬ 
tude  in  a  region  of  sedimentary  cover  is  higher  than  that  in  a 
shield  region.  However,  in  cases  where  the  sedimentary  sections 
contain  limestones  with  high  velocities,  as  at  DH-NY,  the  signal 
can  actually  be  lower  than  in  the  shield  area. 


D .  Analysis  of  TFO  Extended  Array 

During  this  reporting  period  various  studies  have 
been  conducted  to  evaluate  the  use  of  large  arrays,  using  the 
extended  TFO  array  as  the  model.  For  example,  we  have  computed 
the  average  correlations  vs.  distance  and  frequency  for  eight 
eleseismic  signals  and  the  preceding  seismic  noise. 

The  method  of  computing  the  correlations  aligns  the 
P-arrival,  passes  each  trace  through  several  narrow  band  filters, 
and  computes  the  normalized  zero-lag  correlation  coefficients 
between  each  filtered  trace  and  a  standard  trace  at  TFSO.  The 
results  provide  a  measure  of  correlation  vs.  distance  and 
frequency  across  the  extended  array. 

The  average  correlations  vs.  distance  with  the  varia¬ 
tion  in  frequency  as  a  parameter  is  shown  in  Figure  6  for  both 
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Figure  6 


the  IE-SW  line  and  the  NW-SE  line.  These  results  show  that 
the  signal  correlations  are  higher  than  the  noise  correla¬ 
tions  for  almost  all  distances  and  frequencies.  The  noise 
correlations  fall  faster  than  the  signal  correlations  espe¬ 
cially  for  distances  from  3  km  to  100  km.  similar  to  results 
for  single  events  previously  described.  The  correlations  are 
systematically  higher,  the  lowar  the  frequency/  for  both  signals 
and  noise,  as  expected,  since  the  array  dimensions  for  these 
frequencies  are  larger  in  numbers  of  wavelengths. 

All  of  the  signal  correlations  show  a  decided  drop 
between  11  km  and  36  km.  For  distances  less  than  15  km,  the 
correlations  are  between  J-M  seismometer  outputs.  For  distances 
greater  than  15  km,  the  correlations  are  between  J-M  and 
Benioff  seismometer  outputs.  The  two  instruments  differ  in 
both  phase  and  amplitude  response  and  the  difference  is  more 
marked  at  the  higher  frequencies.  It  is  easy  to  show  that 
both  phase  distortion  and  amplitude  distortion  will  cause  a 
lower  correlation  coefficient.  We  suspect  that  the  large  drop 
in  correlations  at  this  range  is  due  to  the  change  of  instru¬ 
ments  and  that  the  signal  correlations  over  the  extended  array 
are  considerably  higher  than  these  results  show.  Since  this 
data  was  analyzed,  arrangements  have  been  made  at  TFO  to  make 
it  possible  to  repeat  the  study  using  all  one  type  instrument. 

We  have  also  made  some  limited  studies  of  travel 
time  anomalies  across  the  extended  array.  The  events  used 
were  18  Kurile  Island  earthquakes. 

All  travel  time  anomalies  were  measured  relative 
to  arrivals  at  MN  NV,  the  nearest  station  in  the  extended  array 
to  the  event  epicenters.  Relative,  rather  than  absolute, 
anomalies  were  measured  since  F-onset  times  for  some  events 
at  some  stations  were  difficult  to  measure.  On  the  other 
hand,  picking  a  later  characteristic  feature  of  the  P-phase 
(a  particular  large  amplitude  peak  or  trough)  was  fairly  easy 
to  measure  anu  provided  a  measure  of  travel  time  intervals 
between  stations. 

The  epicentral  distances  from  MN  NV  ranged  between 
6461  km  and  7564  km.  The  epicenters  were  all  vithin  +_  2.0°  of 
the  same  azimuth. 
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For  the  given  epicenter,  the  Herrin-travel  times 
for  P-arrivals  computed  for  all  stations  were  subtracted 
from  the  observed  arrival  time  of  a  particular  feature  of 
the  P-wave.  Then  these  observed  travel  times,  minus  Herrin 
travel  time  intervals  at  MN  NV,  were  subtracted  from  the 
corresponding  intervals  at  all  other  stations.  The  follow¬ 
ing  Table  1  shows  these  intervals  normalized  to  MN  NV. 

The  travel  time  anomalies  for  IS  events  from  the 
same  epicenter  region  showed  that: 

1.  The  extended  array  responds  fairly  consistent¬ 
ly  to  these  events.  Most  measured  arrival 
times  are  within  +_  C.2  seconds  of  the  average. 

2.  Station  anomalies  are  important.  For  example, 
NL  AZ  is  roughly  0.6  seconds  eaxlier  than 
Kerr in-t ravel  time  predictions. 

3.  Some  of  the  scatter  in  the  arrival  times  for 
these  events  is  probably  due  to  inaccuracies 
in  the  epicenter  locations.  Thus,  the  array 
may  well  respond  more  consistently  than  indi¬ 
cated  by  these  results. 

4.  By  taking  account  of  station  anomalies  and 
expected  travel  times,  the  extended  array 

can  probably  be  velocity-filtered  successfully, 
at  least  for  events  coming  from  the  Kurile 
Islands . 

In  the  future,  this  analysis  will  be  applied  to 
more  events  from  the  Kurile  Islands  and  to  events  from  other 
azimuths . 

E .  Automated  Bulletin  Process 


Under  AFTAC ' s  direction,  an  Automated  Bulletin 
Process  (ABP)  computer  program  was  developed  jointly  by  the 
SDL  and  the  Geotechnical  Corporation.  The  program  associates 
seismic  phase  arrivals  at  a  number  of  stations,  with  epicenters 
reported  by  the  U.  S.  Coast  and  Geodetic  Survey,  and  identifies 
up  to  23  phases.  The  ABP  is  now  being  used  to  accelerate  the 
preparation  of  the  ».jnthly  earthquake  bulletins  of  the  five 
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17  Jul 

12  Aug 

23  Jul 

26  Jul 

27  Jul 

4  Aug 

30  Jun 
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7  Nov 

MNNV 

" 
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+  0.1 
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+  0.1 
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+  0.1 
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+  0.2 
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+  0.1 
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0.0 
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+  0.1 
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0.0 

mm 
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-0.2 

-0.2 

+  0.2 

-0 . 2 
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0.0 
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■EBB 
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+  0.3 
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M 

BO 

+  0.1 

0.0 

0.0 

0.0 

+  0.1 
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-0.7 

-0.3 

-0.3 

-0.5 

-0.4 

-0.7 

-0.6 
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26  Aug 
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0. 

NLAZ 
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-0.1 
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+  C.1 

+  0.2 

0. 
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+  0. 
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+  0.3 
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+  0.4 

+  0.3 
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+  0.8 
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+  0.3 
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+  0.1 
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-0.5 

-0  5 
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-0.2 

-0.5 
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) .  3 
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hfi 
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hi 

hi 

).l 

) .  6 


II  May 


-0,1 

~0.0 

0.0 

0.0 


0.0 
+  0.6 

-0.1 


+  0.2 
-0.2 


VELA-UNIFORM  seismological  observatories  published  by  the 
Geotechnical  Corporation. 

In  terms  of  time  and  manpower  requirements , most 
of  the  work  of  preparing  the  earthquake  bulletins  involves: 

1.  Measurement  of  arrival  times,  amplitude,  and 
period,  preliminary  phase  identification  based 
on  signal  character,  and  card-punching  these 
data  for  each  station.  Analysts  may  record 
trom  300  to  2000  events  per  month  at  each 
station,  depending  upon  the  station,  timr 

of  year,  and  seismic  activity  for  a  given 
month . 

2.  Association  of  phase  arrivals  with  hypocenters 
reported  in  the  U.  S.  Coast  &  Geodetic  Survey 
(USC&GS )  Preliminary  Determination  of  Epicenter 
(PDE)  cards,  identification  of  phases,  merging 
of  associated  and  unassociated  events,  compu¬ 
tation  of  azimuth,  distance,  magnitude,  and 
ground  displacement,  and  finally,  combination 
oE  all  data  from  the  five  observatories  into 

a  form  suitable  for  publication. 

The  ABP  was  designed  to  automate  the  latter  task. 

The  ABP  has  been  written  in  three  parts  (Figure  7). 
These  three  parts  can  be  combined  into  one  program  but  are  usually 
run  independently.  , 

Part  I  computes  the  predicted  arrival  times  of  23 
seismic  phases  from  known  hypocenters  at  each  of  any  number  of 
stations.  Part  II  identifies ,  associates,  and  merges  the  observed 
arrivals  at  each  station  with  the  hypocenters  input  to  Part  I. 

Part  III  rearranges  the  associated  and  unassociated  phases  into 
bulletin  format  and,  by  incorporatinga  program  written  by  J .  G. 
Swanson  of  the  Geotechnical  Corporation,  computes  ground  dis¬ 
placement  and  magnitude.  Part  III  also  computes  distances  for 
local  and  near  regional  events  based  on  tables  of  phase  intervals 
vs.  distance.  Finally,  Part  III  writes  the  entire  bulletin  on 
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Schematic  Routine  On  Which  The  Automatic  Bulletin  Processing  Is  Based 


magnetic  tape  in  the  card  inage  format  from  which  it  is 
published . 


Often,  an  analyst  identifies  a  phase  and  determines 
the  azimuth  and  distance  from  a  hypocenter  to  his  station  from 
the  character  of  the  seismogram  or  by  using  other  criteria. 
Subjective  criteria  are  difficult  to  translate  into  specifications 
for  a  computer  program;  the  analyst  is  encouraged  to  record  this 
data,  even  though  he  has  no  hypocentral  information,  because  in 
the  following  cases  the  analyst's  interpretation  prevails  over 
the  ABP  decision; 

1.  For  events  not  associated  with  hypocenters 
reported  on  the  PD£  cards,  the  analyst's 
identifications  are  used. 

2.  If  the  phase  arrival  is  part  of  a  sequence 
that  has  been  identified  as  a  local  (L) ,  or 
near  regional  (N)  event  by  the  analyst,  and  if 
the  arrival  falls  within  the  expected  arrival 
time  window  of  a  located  teleseismic  (T) ,  or 
regional  (R)  event,  then  the  event  is  left 
unassociated  and  the  station  analyst's  phase 
identifications  are  used. 

3.  The  ABP  does  ,-*ot  attempt  to  identify  Love  or 
Rayleigh  phases.  These  phases  are  identified 
by  the  station  analyst. 

In  all  other  cases,  the  ABP  overrides  or  blanks  out 
the  analyst's  identification.  However,  the  choice  of  which  of  the 
analyst's  phase  identifications  the  program  will  supercede,  can 
be  specified. 


Allowable  residuals  for  phase  identification  are 
based  on  statistical  distribution  of  time  residuals  observed  for 
each  of  the  23  phases  which  ABP  attempts  tc  identify.  A  digital 
magnetic  tape  was  prepared  containing  focal  depth,  station-to- 
epicenter  distance,  station  name,  and  phase  arrival  time  for 
all  identified  phases  from  1613  known  hypocenters  recorded  at 
BMSO,  CPSO,  TFSO,  UBSO ,  and  WMSO  during  April,  May,  June,  and 
October  1963 . 


From  these  data,  the  23  most  commonly  recorded  phases 
were  selected,  the  mean  arrival  time  residuals  were  determined 
for  those  phases,  and  a  "maximum  allowable"  residual  value,  was 
established  for  each  phase,  based  on  commonly  observed  frequencies 
of  the  phase.  The  reamining  phases  such  as  PcS ,  PSS ,  SSP ,  etc., 
represe,  :  less  than  one  per  cent  of  the  total  number  of  phases 
identified  by  the  analysts. 

Any  residuals  which  fell  outside  this  "maximum 
allowable"  value  were  assumed  to  be  the  result  of  inaccuracy 
of  the  PDE  data,  data  transcription  errors,  or  questionable 
phase  identifications,  and  were  not  used  in  the  determination  of 
the  mean  residual  for  that  phase.  All  association  windows  are 
approximately  equivalent  to  the  mean  arrival  time  residual  plus 
or  minus  twice  the  standard  deviation,  except  in  the  cases  of 
P  and  PKP  where  the  windows  are  equivalent  to  the  mean  residuals 
plus  or  minus  three  times  the  standard  deviation. 

The  ABP  takes  six  hours  on  the  CDC  1604  co  process 
one  month's  data  for  the  five  VELA  observatories.  Each  station's 
arrivals  are  processed  independently  and  merged  with  data  already 
processed  for  other  stations.  Under  this  system,  a  bulletin 
can  be  generated  for  any  number  of  stations,  or  the  program  can 
be  restarted  after  a  one,  two,  or  three  station  bulletin  has 
been  written  on  magnetic  tape  for  a  given  month. 

To  date,  the  ABP  has  been  run  for  the  September, 
October,  November,  and  December  1964  bulletins. 

Tulletin  data  from  February  ’963  through  December 
1964  has  been  written  on  magnetic  tape  in  a  compact  format  to 
facilitate  retrieval.  A  set  of  programs  has  been  written  to 
process  the  Bulletin  data,  summarizing  the  frequency  of  occurrance 
of  each  phase  at  each  observatory  as  a  function  of  earthquake 
magnitude,  depth,  and  distance. 

F .  Shot  and  Earthquake  Analyses 
1 .  SALMON 

SALMON,  the  first  shot  of  the  DRIBBLE  series,  was 
detonated  at  the  Tatum  Salt  Dome  near  Hattiesburg,  Mississippi, 
on  22  October  1964.  Thirty-eight  stations  recorded  short- 
period  signals  and  20  stations  recorded  long-period  signals. 


20 


A  Mexican  earthquake,  occurring  approximately  4  minutes  after 
the  SALMON  explosion,  possibly  interfered  with  western  U.  S. 
stations  recording  the  surface  waves  from  SALMON.  The  average 
magnitude  for  SALMON  was  4.58,  and  tv^.  stations,  Eutaw,  Alabama 
and  Jena,  Louisiana,  showed  compressional  first  motion  as 
defined  by  the  TWG  II  First  Motion  Criteria.  The  travel-time 
residuals  from  the  Pn  and  P  phase  corresponded  favorably  with 
those  derived  from  the  GNOME  event  and  showed  that,  generally, 
the  apparent  upper  mantle  velocities  are  substantially  higher 
to  the  east  and  north  than  to  the  west  .owards  the  Nevada  Ter;t 
Site.  The  most  distant  station  that  recorded  SALMON  was  La  ^az, 
Bolivia,  at  a  distance  of  5704  kilometers. 

2 .  Texas-Louisiana  Earthquake  (24  April  1964) 

The  average  magnitude  was  4.36,  and  short -period 
signals  were  recorded  by  30  stations;  no  long-period  signals  were 
recorded  due  to  interference  from  a  larger  magnitude  (6.3)  New 
Guinea  earthquake.  Two  stations  Jena,  Louisiana  and  Grapevine, 
Texas,  showed  rnrefactional  first  motion  as  defined  by  the 
TWG  II  First  Motion  Criteria. 

3 .  West  Virginia  Earthquake  (25  November  1964) 

The  average  magnitude  was  3.95,  and  short -period 
signals  were  recorded  by  27  stations;  the  station  at  Beckley, 

West  Virginia  recorded  a  long-period  surface  wave,  identified 
as  a  high  frequency  sedimentary  Rayleigh  wave.  None  of  the 
stations  showed  compressional  or  raref act ional  first  motion 
as  defined  by  TWG  II  First  Motion  Criteria. 

4 .  Seismic  Waves  From  the  SS  Village  Explosion 

This  was  a  chemical  explosion  at  approximately 
7,000  feet  underwater  wh'ch  took  place  aboard  a  ship  being 
sunk  in  the  Atlantic  Ocean  off  the  coast  of  New  Jersey  and 
Delaware  on  17  September  1964.  The  average  magnitude  was  4.45, 
and  short -period  signals  were  recorded  by  30  stations;  no 
stations  recorded  any  long -period  data.  Four  stations,  Delhi, 

New  York,  Beckley,  West  Virginia,  Houlton,  Maine,  and  Jerome, 
Arizona,  showed  compressional  first  motion  and  one  station, 
Berlin,  Pennsylvania,  showed  rarefact ional  first  motion,  as 
defined  by  the  TWGII  First  Motion  Criteria. 
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The  character  of  the  records  was  definitely 
influenced  by  r  ltiple  water  reflections  which  appeared  to 
repeat  a  numbei  of  times.  Although  the  time  intervals  between 
peak  amplitudes  appeared  about  constant,  the  onsets  of  the 
individual  reflections  were  difficult  to  pick.  It  was  also 
noted  that  the  character  of  the  multiples  varied  from  station 
to  station  indicating  possible  receiver  and  path  effects. 

In  Figure  8,  three  typical  seismograms  are 
shown;  time  intervals  of  the  Pn  and  P  multiples  as  determined 
are  consistant  with  those  from  many  other  records.  These  are 
indicated  by  the  solid  lines.  The  spacing  of  the  multiples 
(2.6  seconds)  is  dependent  on  the  thickness  of  the  water 
layer.  With  the  hypothetical  shot  position  as  shown,  the 
shot  should  produce  another  set  of  pP  multiples  equal  in 
amplitudes  to  the  corresponding  P  multiple  amplitudes.  These 
are  not  readily  apparent;  so  it  is  probable  the  shot  occurred 
at  the  liquid-solid  interface.  However,  the  absence  of  a 
well  defined  pP  multiple  could  also  be  explained  by  distor¬ 
tion  of  the  wave  form  caused  by  the  time  lag  and  phase  re¬ 
lationships  of  the  multiple  with  some  other  phase. 

G .  REMODE  (Rectili tear  Motion  Detection) 


Studies  *  rectilinear  motion  detection  at  SDL 
evolved  from  a  long  series  of  experiments  related  to  the  use  of 
the  analog  computer  in  processing  seismograms  recorded  on 
magnetic  tape  in  analog  form.  These  studies  have  led  to  a 
class  of  operators  which  we  have  named  REMODE,  standi »ig  for 
Rectilinear  Motion  Detection.  Our  REMODE  repertoire  contains 
both  analog  and  digital  programs,  both  of  which  we  are  actively 
testing  and  modifying  as  better  techniques  are  discovered. 

Early  studies  of  this  series  included  the  sum  of 
the  squares  of  three-component  seismograms  on  which  we  have 
already  reported.  Later  studies  included  matched  filters. 

A  signal  waveform  which  is  reasonably  constant 
from  place  to  place  for  a  given  type  of  event  such  as  an 
explosion  would  make  possible  the  use  of  a  matched  filter 
to  enhance  signal  detection  in  an  automatic  system.  The 
purpose  of  this  study  was  two-fold: 
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Figure  8 


1.  To  determine  to  what  extent  signal  waveforms 
from  explosions  are  consistent;  and, 

2.  To  determine  whether  a  matched  filter  success¬ 
fully  enhances  the  determination  of  first  motion. 

Similarities  of  waveforms  have  been  determined  for 
two  Russian  and  two  Algerian  events  as  recorded  at  BL  WV,  DR  CO, 
and  LC  NM.  Th^  waveforms  were  presented  as  Laguerre  approxi¬ 
mations  (20  coefficients)  to  the  initial  pulse.  Averages  for 
the  signals  from  each  region  were  also  made.  Comparison  of  the 
waveforms  s howed  considerable  similarity  but  also  some  dis¬ 
advantages  of  this  method  of  approximation.  For  example,  many 
of  the  individual  waveforms  displayed  a  high  frequency  noise 
introduced  by  the  approximation. 

The  average  signal  for  two  evnets  analyzed  was 
used  as  a  matched  filter  to  process  seismograms  from  three 
events  as  recorded  at  BL  WV,  DR  CO,  and  LC  NM.  In  each  case, 
the  third  seismogram  was  noisier  than  the  two  for  which  the 
signal  approximations  were  derived. 

The  original  data  traces  were  processed  through 
the  matched  filter  formed  from  the  average  signal.  The  data 
traces  were  then  processed  through  the  quadrature  filter 
(same  amplitude  response  as  the  matched  filter  but  with  a 
90°  phase  shift).  If  we  square  each  of  the  corresponding 
filtered  traces  and  add,  the  resulting  correlation  traces 
display  sharp  maxima  related  to  the  original  signals  as  shown 
in  Figure  9.  The  time  origin  cnoeen  for  the  average  waveform 
is  such  that  the  sharp  maxima  on  the  final  traces  indicate 
the  onset  of  the  signal.  The  following  table  shows  the 
precision  found  for  these  events  and  stations. 
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Results  Of  Squaring  And  Summing  The  Direct  And  Quadratured  Matched  Filter 
Outputs 


Figure  9 


TRAVEL  TIMES  DETERMINED  BY  MATCHED  FILTER  TECHNIQUE 

FOR  SAHARAN  EVENTS 


DATE 

BL  WV 

DR  CO 

LC  NM 

1  May  1962 

Detector  Time 
Origin  Time 

10:11:31.2 
10: 00: 00.3 

10:13:12  .2 
10:  00:00.3 

10:13:21.0 

10:00:00.3 

Travel  Time 

11:30.9 

13:11.9 

13:20.7 

18  March  1963 

Detector  Time 
Origin  Time 

10:13:31.2 

10:02:00.5 

10:15:12  .4 
10:02:00.5 

10:15:21.4 
10:  02: 00.5 

Travel  Time 

11:30.7 

13:11.9 

13:20.9 

20  October  1963 

Detector  Time 
Origin  Time 

13:11:30.9 

13:00:00.1 

13:13:12.0 

13:00:00.1 

13:13:20.8 

13:00:00.1 

Travel  Time 

11:30.8 

13:11.9 

13:20.7 

In  the  study  of  matched  filters  it  was  learned  that 
there  was  a  danger  of  false  alarms  which  might  be  overcome  or 
reduced  by  some  sort  of  filtering.  Frequency  filtering  will  not 
satisfy  this  requirement  because  the  false  alarms  generally  have 
their  origins  in  the  frequency  band  of  the  matched  filter.  As 
a  solution,  we  have  considered  coherency  filtering  in  several 
forms.  Further,  in  order  to  better  understand  the  processes 
being  studied,  we  have  examined  coherency  filtering  of  very 
narrow  passbands . 

We  first  made  up  pairs  of  synthetic  seismograms 
by  mixing  a  given  signal  with  different  samples  of  noise. 

Each  seismogram  trace  is  decomposed  into  seven  narrow- 
band  component  s : 

Z  (t)  =  Z(t)*W  (t)  and  R  (t)  =  ?(t)*W  (t) 
n  n  n  n 

where  n  varies  from  one  to  seven.  The  filters  in  this  experi¬ 
ment  were  constructed  with  Q  =  5,  and  the  center  frequencies 
were  0.64,  0.8,  1,  1.25,  1.6,  and  2  cps. 
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For  any  given  frequency  component,  the  simulated 
vertical  component  was  correlated  with  the  simulated  horizontal 
component  to  get  the  correlation  function: 


C  (t) 
n 


Z  (t)  *  R  (t) 

n  n 


rz  2 

l  n 


(t) 


*5 


"Cn(t) 

A  gain  function  G  (t)  of  the  form  e  was  constructed  and 

n 

made  to  operate  on  the  vertical  trace.  The  final  output  was 
composed  by  summing  the  seven  components 


7 

Z* (t )  =  V  G  (t)  Z  (t) 
n  n 

n=l 


Similar  correlations  were  made  basing  the  correlation 
criteria  and  gain  control  G'  on  the  comparative  amplitudes  one 
would  expect  to  get  on  the  two  components  for  a  given  angle  of 
emergence  or  azimuth.  We  also  performed  cross-correlation  between 
traces  from  adjacent  frequency  bands  and  used  the  resulting 
function  as  a  gain  control  on  one  of  the  traces.  The  actual 
cross-correlation  function  usr  ... .  given  frequency  band  was 
the  average  of  the  correlations  between  the  trace  for  that  band 
and  the  traces  for  adjacent  bands  on  either  side  (where  both 
exist ) : 


Z*(t)  *  Z*  (t)  ZMt )  *  Z*  (t) 

n _ n-1  +  n _ n+ 1 _ 

Zn2(t)  ‘  Z?,-l  k2(t)  '  Z*n+1  (t)  3  "2 

-c  (1 

and  the  corresponding  gan  control  G^(t)  was  of  the  form  e  n 

An  example  of  the  output  for  a  single  frequency  band 
is  shown  in  column  B  of  Figure  10;  the  band  shown  is  the  center 
band  of  the  seven  used  and  is  that  of  the  apparent  frequency 
of  the  input  signal.  In  column  C  of  the  same  figure  are  shown 
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examples  of  the  final  output  defined  by 
7 

2;***  (t)  =Y  G"  ( t )  G  ( t )  G  '  ( t )  Z  (t) 
L  n  n  n  n 

n=l 


We  now  consider  the  basis  of  REMODE  to  lie  in  the  fact 
that  the  onset  of  P-wave  signals  tends  to  be  in  phase  on  both  the 
horizontal  and  vertical  components  of  seismograms;  conversely, 
seismic  noise  tends  to  be  randomly  in  or  out  of  phase  on  the 
horizontal  and  vertical  components.  This  relationship,  as  well 
as  others,  is  useful  as  a  basis  for  extraction  of  seismic 
signals  from  noise.  Specifically,  for  example,  a  raypath 
corresponding  to  a  body  wave  phase  initially  emerges  at  the 
receiver  at  a  unique  angle  depending  largely  on  the  distance 
from  the  source.  Such  an  arrival  is  expected  to  produce 
similar  particle  motion  on  the  radial  and  vertical  components. 

On  the  other  hand,  the  horizontal  and  vertical  components  of 
seismic  noise  are  often  observed  to  have  very  little  simi¬ 
larity  in  waveform;  the  noise  field,  consisting  of  many 
vibrational  modes  of  different  orbit  and  direction  of  propa¬ 
gation,  can  be  expected  to  produce  rapidly  changing  and  largely 
uncorrelated  waveforms  on  both  the  horizontal  and  vertical 
components . 

A  method  is  drawn  from  electrical  engineering  to 
determine  the  operator  which  extracts  a  least  squares  estimate 
of  a  common  waveform  present  in  two  time  series.  For  weak 
body  waves,  the  application  of  the  extraction  operator  improves 
the  signal-to-noise  ratio,  and  isolates  phases  which  cannot 
otherwise  be  detected  with  a  narrow  band-pass  filter.  Strong 
body  wave  phases,  which  would  otherwise  be  obscured  by  signal 
generated  noise,  are  also  isolated. 

We  distinguish  between  two  classes  of  operators: 
zero  lag  correlation  coefficients  are  labelled  "scalar  para¬ 
meters"  and  correlation  functions  involving  the  convolution 
integral  are  labelled  "functional  parameters".  In  general, 
where  as  prior  work  in  this  fiel  1  has  been  restricted  to  use 
of  the  zero  lag  correlation  coefficients,  we  recommend  use  of 
the  more  comprehensive  correlation  functions.  In  particular, 
we  have  found  the  best  results  in  using  as  a  gain  control  the 
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even  part  of  the  cross-correlation  function  between  the  two 
traces  in  question. 

Figure  11  illustrates  the  efficiency  of  the  even  part 
of  the  correlation  function  in  extracting  signal  from  Gaussian 
noise . 


We  have  mixed  a  given  signal  with  different  random 
noise  sampler  at  a  peak  signal-to-noise  ratio  of  one.  Gaussian 
noise  was  obtained  from  a  random  noise  generator  passed  through 
a  linear  filter  which  simulates  the  Benioff  system  response. 

The  processed  seismograms  are  shown.  In  each  case, 
the  operator  r amber  indicates  the  trace  on  which  the  filtering 
was  performed.  It  is  evident  that  a  considerable  improvement 
has  been  attained  in  the  signal-to-noise  rrcio,  at  the  same 
time  maintaining  the  signal  shape  reasonably  .ell. 

By  contrast,  the  results  of  band-pass  filtering  on  the 
same  synthetic  records  are  distinctly  inferior  and  one  would 
not  expect  to  improve  even  after  considerable  effort  to  find 
the  optimum  band-pass  filter. 

We  feel  that  small  arrays  of  three-component  seis¬ 
mometers  are  the  most  suitable  for  REMODE  filtering.  As  a 
preliminary  trial  we  present  an  example  of  REMODE  filtering  on 
the  7-element  array  at  Oslo,  Norway  (Figure  12).  It  can  be 
seen  that  an  almost  unpickable  seismogram  has  been  converted 
to  a  useful  record  by  the  REMODE  operation.  In  this  case,  since 
there  was  only  one  three-component  instrument  at  the  site,  the  F 
vertical  seismograms  were  velocity  filtered  and  summed  before 
REMODE  filtering. 

The  actual  shape  of  the  pulse  passed  by  the  filter 
is  likely  to  be  quite  distorted,  especially  when  the  S/N  ratio 
is  low,  since,  in  addition  to  distortions  introduced  by  the 
filter  itself,  the  sum  of  the  noise  and  signal  is  admitted, 
inasmuch  as  the  filter  acts  somewhat  like  a  gate.  Hence, 
this  filter  provides  information  on  the  time  of  first  motion 
and  the  interval  of  rect ilinearity  for  a  phase  but  little 
highly  reliable  information  on  the  pulse  shape.  Indeed  one 
might  expect  the  signal  as  admitted  by  a  REMODE  filter  to 


27 


SE«SMOGRAMS 


Z  l 

z  z 

21 

l  4 

Z5 

Z6 

Z  7 

-AIM  2’s 

T 

REMQOE  PASS) 

F‘LTrRrO 
S-V  1 

F<_TL»ro  T 

RE  MODE  (2nd  PASS) 

ftlTEPED 
"UV  4 

;*LTFwt!D  T 

time 


u*  A 

7  \r:  f— 


'AJ\, 


■v 


‘r  \v/-v-v'- 


pl  <»  PcPIJB) 


T  * 

♦ 

16 II 00  oz 


REMODE  Array  Processed  oeisroograms 


for  OSLO,  NORWAY, 


a  =  7565  Km, 


SALMON  Event 


Figure  12 


look  much  like  a  section  of  the  noise  when  the  S/N  ratio  is 
low.  This  appears  to  be  the  case  with  the  REMODE  trace  for 
the  Oslo  array,  although  the  true  shape  of  the  pulse  is  itself 
unknown . 


H .  Supp1 emental  Analysis  of  SALMON  Data 

A  rectilinear  motion  detection  (REMODE)  technique  has 
been  used  to  isolate  arrivals  from  the  SALMON  event  at  tele- 
seismic  distances  in  several  cases  where  the  signal-to-noise 
ratio  was  low  and  visual  methods  or  the  ordinary  array  pro¬ 
cessing  procedures  failed.  The  essential  aspects  of  array 
and  single  station  REMODE  processors  currently  in  use  were 
outlined  and  examples  from  the  SALMON  event  were  processed. 

The  phase  isolation  capability  cf  the  processor  seemed  most 
promising  and  was  fully  exploited  in  the  analysis  of  the 
SALMON  data.  In  particular,  P-wave  travel  time  data  was  ob¬ 
tained  from  the  processed  traces  along  three  azimuths  from 
the  SALMON  site.  The  travel  time  curves  obtained  showed 
reasonably  well-defined  multiple  branches.  An  ir/ersion 
method  giving  the  velocity  structure  appropriate  to  these 
data  was  briefly  described,  and  the  structures  obtained 
from  the  SALMON  data  were  discussed.  Radiation  patterns 
appropriate  to  the  primary  travel  time  branch  for  the  P-wave 
were  obtained  for  various  frequencies,  and  these  were  interpreted 
in  terms  of  an  ideally  symmetric  and  purely  compressional  source. 
Comparisons  were  also  made  to  theoretically  derived  radiation 
patterns  for  explosive  sources  with  associated  tectonic  energy 
release . 


The  surface  wave  dispersion  along  the  three  profiles 
from  the  SALMON  site  has  been  measured  and  an  inversion  method 
for  this  data  was  discussed.  Velocity  structures  appropriate 
to  the  surface  wave  data  were  described  and  compared  to  the 
body  wave  structures.  Rayleigh  wave  radiation  patterns  were 
also  obtained  and  interpreted  in  terms  of  a  compressional 
source . 


A  few  general  observations  concerning  the  source  and 
the  associated  seismic  field  can  be  made  from  a  preliminary 
analysis  of  the  signal  from  this  event. 

Fundamental  mode  Rayleigh  waves  were  identified  at 
21  LRSM  stations,  in  the  ^cr?od  range  8  to  40  sec.,  out  to 
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about  2500  km.  Higher  mode  Rayleigh  waves  were  observed  to 
distances  of  3000  km  on  the  continent  and  at  nearly  every 
operative  LRSM  station  P  phases  were  detected  at  all  but 
one  or  two  continental  LIv  stations  and  at  La  Paz,  Bolivia, 
at  a  distance  of  5704  km.  The  long  period  Rayleigh  wave 
excitation  was  not  very  great,  as  might  be  expected,  and  the 
signal-to-noise  ratio  was  generally  quite  low,  except  at  the 
closest  pair  of  staions .  In  addition  the  occurrence  of  an 
earthquake  in  Mexico  either  complicated  or  completely  obliterated 
the  Rayleigh  waves  at  some  of  the  western  stations.  Consequently, 
the  analysis  of  surface  wave  radiation  was  limited  to  14  stations. 
On  the  other  hand,  the  quality  of  the  short  period  data  was 
uniformly  good  and  uncomplicated  by  other  events,  so  that  all 
LRSM  stations  (40)  were  used  in  the  body  wave  analysis. 

As  with  the  GNOME  event,  no  Love  wave  radiation 
was  observed.  This  is  in  contrast  to  most  of  the  NTS  shots, 
for  example,  SHOAL  and  BILBY,  which  produce  SH  surface  waves 
of  this  type.  One  may  note  that  the  presence  of  Love  waves 
would  be  expected  if  tectonic  energy  release  occurred.  While 
NTS  is  in  a  tectonically  active  area,  and  the  medium  would 
normally  be  prestressed,  SALMON  and  GNOME  were  in  salt  and 
would  not  support  even  a  modest  prestress  field.  Thus,  the 
lack  of  SH  radiation  from  SALMON  and  GNOME  would  be  consistent 
with  the  presence  of  such  radiation  from  NTS  events  if  tectonic 
release  were  involved.  The  observations  are  still,  however, 
open  to  other  interpretations. 

P-wave  arrivals  along  the  NW  profile  were  approxi¬ 
mately  the  same  as  the  J-B  times,  while  on  the  N  profile  they 
were  4  to  5  seconds  earlier  than  the  J-B  times.  These  ob¬ 
servations  were  in  general  agreement  with  those  from  GNOME. 

Figure  13  shows  two  especially  good  examples  of  the 
signal-to-noise  enhancement  and  phase  isolation  capability  of 
the  filter.  The  traces  shown  correspond  to  single  LRSM  stations, 
so  chat  the  array  processor  degenerates  to  a  simple  double 
application  of  the  REMODE  filter  in  this  case.  The  event  shown 
is  SALMON,  and  the  stations  are  in  the  20°  distance  range  and, 
as  might  be  expected,  several  clear  P  arrivals  have  been  isolated, 
presumably  corresponding  to  an  inverse  branch  or  multiple  branches 
of  the  travel  time  curve.  As  was  previously  pointed  out,  these 
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REMODE  Seismograms  for  SALMON  (EK-NV  and  RY- 


filtered  traces  should  correspond,  for  a  properly  designed 
REMODE  filter,  to  the  nearly  rectilinear  first  motions  of  the 
individual  phases. 

Of  the  stations  processed  in  this  study,  roughly 
30  percent  gave  REMODE  outputs  with  S/N  enhancement  and  phase 
isolation  comparable  to  the  examples  shown. 

The  most  common  failing  in  the  remaining  processed 
records,  aside  from  occasional  reductions  in  the  overall  signal 
to  noise  enhancement  feature  of  the  process,  was  the  absence 
of  a  clearly  isolated  Pn  arrival  in  the  zone  where  the  S/N 
ratio  for  Pn  is  low.  While  this  apparent  lack  of  sensitivity 
may  be  due  to  the  manner  of  data  presentation  (i.e.  the  scaling 
of  the  plotted  filter  output),  it  most  likely  indicates  a  threshold 
of  S/N  for  detect ibility ,  with  this  particular  filter  at  least. 

The  precise  thre  hold  limits  are  observed  to  vary  from  .location  to 
location  however,  so  that  one  must  conclude  that  the  character 
of  the  noise  and  signal,  but  especially  of  the  noise,  determines 
to  a  great  extent  the  sensitivity  of  the  filter  to  small  ampli¬ 
tude  rectilinear  signals. 

In  examining  signal  spectra  as  a  function  of  distance, 
it  is  clear  that  the  maximum  spectral  amplitude  for  P  shifts  to 
longer  periods  with  increasing  distance,  due  to  the  roughly 
linear  increase  in  attenuation  with  frequency.  The  observed 
attenuation  will,  of  course,  be  affected  by  thg  variation  of  the 
dissipation  function  Q  with  depth  in  the  earth  so  that  the 
attenuation  will  be  dependent  on  the  actual  minimum  time  path 
followed  by  the  phase.  This  effect  would  be  of  second  order, 
however,  compared  to  the  much  larger  amplitude  affects  associated 
with  the  velocity  variations  with  depth,  particularly  with  velocity 
reversals  and  rapid  or  discontinuous  changes  in  velocity. 

Thus,  the  observed  spectra  would  not  be  expected  to 
fall  off  uniformly  with  distance  in  the  range  where  the  minimum 
time  paths  would  be  strongly  affected  by  the  velocity  variations 
in  the  upper  mantle  of  the  earth.  Indeed,  strong  and  rapid 
variations  with  distance  would  be  expected  in  view  of  the 
likelihood  of  large  gradients  in  the  velocity  variation  with 
depth  and  the  existence  of  the  low  velocity  zone  in  the 
upper  mantle. 
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A  comparison  of  spectral  amplitude  with  distance 
at  particular  frequencies  is  given  in  Figure  14.  The  anomalous 
behavior  only  suggested  by  the  individual  spectra  is  now  quite 
apparent.  The  essential  aspects  of  the  amplitude-distance 
variation  obtained  and  illustrated  in  this  figure  have  been 
observed  many  times  before  by  otner  investigators,  and  in 
particular,  for  events  observed  along  a  profile  of  stations 
in  this  same  region  .  The  spectral  detail  in  the  present  case 
is  perhaps  somewhat  greater  than  has  been  available  before, 
so  that  a  correspondingly  more  detailed  interpret  it  ion  will  be 
advanced . 


The  spectra  will  be  interpreted  together  with  the 
observed  travel  time  curves.  However,  viewed  independently, 
they  show  a  fairly  rapid  fall-off  with  distance  out  to  approxi¬ 
mately  1250  kilometers.  The  rather  gradual  increase  in  ampli¬ 
tude,  beginning  at  this  point  and  reaching  an  apparent  maximum 
at  1800  to  1900  kilometers,  probably  indicates  a  first  arrival 
P-phase  different  from  Pn .  A  decrease  in  amplitude  at  around 
2000  km,  followed  by  a  second  peak  in  the  spectrum  again  probably 
indicates  a  second  P-phase. 

This  rather  superficial  interpretation  has  been  checked 
with  the  travel  time  observations  for  evidence  of  changes  in 
the  first  arrival  curve  slope  as  well  as  for  the  existence  of 
multiplicities  or  cusps  in  the  curves  accounting  for  the 
amplitude  maxima. 

Figure  15  shows  the  travel  time  versus  distance 
data  for  all  profiles  combined.  The  d~ta  shown  correspond  to 
first  arrival  phases  and  all  phases  within  20  seconds  of  the 
first  arrival  at  a  given  distance,  with  the  exception  of  those 
identified  as  PP  and  PPP .  At  the  greater  distances  some  of  the 
later  arrivals  are  shown  since  they  were  often  of  large  amplitude 
and  may  be  related  to  mantle  velocity  variations  which  manifest 
themselves  in  the  first  arrival  curves  at  somewhat  greater  dis¬ 
tances  than  those  considered  in  the  present  study.  The  travel 
time  observations  obtained  for  PP  and  PPP  from  the  REMODE  pro¬ 
cessed  time  series  have  the  expected  character  as  functions 
of  distance  and  serve,  in  part,  to  increase  the  confidence 
with  which  one  may  interpret  the  travel  times  for  other  phases, 
particularly  the  P  and  P  travel  times. 

The  solid  circles  in  the  figures  indicate  large  amplitudes 
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arriving  at  a  particular  distance,  open  circles  relatively  small 
or  moderate  amplitudes.  The  other  symbols  refer  to  J-B  travel 
times . 


The  travel  time  curves  drawn  through  the  data  are  based 
on  the  travel  time  data  themselves  and  the  previously  discussed 
amplitude-distance  spectra.  In  addition,  the  curves  correspond 
in  their  general  characteristics  to  theoretical  curves  obtained 
from  the  velocity  structure  shown  in  Figure  16.  The  initial 
selection  of  this  structure  was  based  upon  a  consideration  of 
the  long  period  surface  wave  dispersion  data  obtained  by  a  ^ 
number  of  investigators  and  compiled  by  Anderson  and  Tokstiz. 

The  velocity  structure  in  the  mantle  is  a  modified  version  of 
the  CIT  11  oceanic  model  proposed  by  Anderson  and  Toksttz  to 
explain  the  Love  wave  dispersion  data.  The  model  has  also  been 
found  to  provide  good  amplitude  and  travel  time  agreement  with 
independent  body  wave  data. 

The  only  extensive  modification  to  their  structure  is 
the  addition  of  a  crustal  section  -  and  of  course,  conversion  to 
compressional  velocities,  which  involves  assumptions  of  Poisson's 
ratios.  Thus,  in  view  of  the  essential  agreement  between 
continental  and  oceanic  dispersion  for  periods  beyond  about 
50  seconds,  the  model  agrees  with  continental  long  period 
surface  wave  dispersion. 

While  the  structure  has  been  checked  and  perturbed 
somewhat  for  consistency  with  the  P  phases  observed,  no 
particular  effort  has  been  made  to  adjust  the  fit  to  the 
Pn  and  other  Dody  phase  data  except  in  an  approximate  manner. 

For  this  reason,  the  structure  shown  should  be  regarded  as 
a  preliminary  fit,  which  must  be  adjusted  to  account  more 
fully  for  Pn,  the  reflected  P  phases,  and  the  finer  details 
of  the  direct  P  phases.  In  addition,  the  crustal  section 
used  here  corresponds  to  an  average  crust  which,  for  a  finer 
fit,  would  have  to  be  replaced  by  the  crustal  section  appropri¬ 
ate  to  the  particular  region  or  regions  involved  in  the  travel 
time  profile. 
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WORK  IN  PROGRESS 


Work  is  continuing  on  several  of  the  projects  discussed 
in  Section  II.  In  addition  to  these  projects,  other  work  in 
progress  at  the  SDL  includes  the  following  items: 

A .  Inversion  of  Surface  Wave  Phase  and  Group  Velocity 

Dispersion  Observations 

Preliminary  results  using  the  comb-filter  program 
fo  'asuring  surface  wave  group  velocity  indicate  that  this  is 
a  !•'  >rful  and  sensitive  tool  for  observing  crustal  and  upper 
ma  ^  i  structure.  In  order  to  utilize  fully  the  potential  of 
the  method,  an  analytical  procedure  for  inverting  the  obser¬ 
vations  to  find  the  corresponding  best-fitting  structure  is 
required . 

Following  a  suggestion  of  Jeffreys  (1961)^,  Archambeau 
and  Andersen  ^  derived  the  inversion  coefficients  for  Love  wave 
phase  velocity  observations,  in  terms  of  the  potential  and 
kinetic  energy  content  of  the  layers  comprising  the  structural 
mode.  Flinn  has  recently  derived  inversion  coefficient  for 
Rayleigh  wave  phase  velocity,  again  in  terms  of  layer  energies. 

If  the  observations  are  taken,  holding  frequency  instead  of  wave 
number  constant  (as  the  comb-filter  group  velocity  approach 
does),  the  inversion  coefficients  are  layer  energy  fluxes. 

As  far  as  phase  velocity  is  concerned,  the  Rayleigh  wave 
inversion  coefficients  are  exactly  as  one  expects  physically: 
the  layer  potential  energy  is  split  into  two  parts,  a  com- 
pressional  part  which  affects  only  the  variation  of  layer  bulk 
modulus,  the  other  part  rotational,  which  affects  only  the 
layer  rigidity. 

Flinn  and  Archambeau  derived  expressions  for  the 
inversion  of  Love  wave  group  velocity  observations,  again  in 
term  oi  layer  energies  or  energy  fluxes.  The  inversion  coef¬ 
ficients  for  Rayleigh  wave  group  velocity  are  much  more  com¬ 
plicated,  however,  and  no  expression  for  the  corresponding 
coefficients  in  terms  of  layer  energies  has  yet  been  found. 
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B .  Perturbation  Theory  for  the  Inversion  of  Body 
Wave  Travel  Time  Data 


Estimation  of  the  internal  elastic  structure  of  the 
earth  from  body  wave  travel  times  ordinarily  uses  both  arrival 
time  data  and  the  slope  of  a  curve  fitted  to  the  observed  travel 
times.  There  are  well  known  difficulties  in  fitting  such  curves, 
particularly  in  cases  where  the  observed  travel  times  suggest  a 
multivalued  travel  time  function.  A  perturbation  method  has 
been  developed  in  which  an  initial  velocity  distribution,  which 
may  be  chosen  to  be  in  agreement  with  surface  wave  dispersion 
data  over  the  region  in  question,  is  pertubed  iteratively 
unti.l  the  theoretical  travel  time  function  is  in  agreement,  in 
the  least-squares  sense,  with  the  observed  data.  The  method 
does  not  require  any  estimate  of  the  curve  slope  or  multi¬ 
plicity,  nor  does  it  require  a  particularly  dense  and  complete 
data  coverage  over  the  entire  distance  range.  Several  phases 
can  be  used  simultaneously. 

The  initial  velocity  distribution  is  taken  to  be  in 
the  form  V(r)  =  Ea  f  (r) ,  where  the  f  are  any  set  of  finite 
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differentiable  functions  over  the  appropriate  depth  range.  At 
each  step  in  the  iteration  process,  variations  in  the  are 

computed  from  the  travel  time  residuals  6T  at  fixed  distance, 
from  the  relation 


where  T  represents  a  contribution  to  travel  time  due  to  the 
term  a^f^(r) . 

C .  Recursive  Numerical  Filters 

The  objective  of  this  study  is  to  design  a  digital 
computer  program  that  will  taxe  conventional  filter  parameters 
and  automatically  design  an  equivalent  recursive  digital  filter. 
The  parameters  used  are  the  poles  and  nume  itor  polonomial 
coefficients  in  the  transfer  function  of  the  conventional 
filter.  The  resultant  digital  filters  are  exact  in  the  sense 
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that  their  impulse  responses  are  equal  to  the  uniform  samples 
of  the  impulse  responses  of  their  prototype  analog  filters. 

More  concisely,  the  objective  is  to  write  a  computer  program 
to  convert  Laplace  transforms  to  Z-transforms . 

Several  intermediate  and  related  results  have  been 
obtained.  These  include  the  following: 

1 .  A  program  to  compute  the  frequency  response  of 
a  conventional  filter  has  been  written.  This 
program  is  used  as  a  check  on  the  digital  model 
in  the  frequency  domain  and  has  been  used  to 
check  theoretical  versus  measured  frequency 
responses  for  seismometer  systems. 

2.  A  program  to  design  quadrature  filters  has  been 
written.  These  filters  are  needed  to  compute 
analog  (and  eventually  digital)  envelope  functions. 
Results  using  these  filters  on  the  analog  computer 
are  very  good . 

3.  Recursive  models  of  the  elements  in  a  seismic 
system  have  been  developed. 

4.  Various  recursive  Chebychev  and  Butterworth 
low-pass  filters  were  designed  and  compared  with 
their  prototypes.  Results  were  good,  in  general, 
and  an  attempt  is  underway  to  improve  those 
cases  which  were  not  satisfactory. 

D .  Correlogram  Analysis  from  Linear  Arrays 

In  the  past  six  months  correlograms  from  over  a  hundred 
teieseismic  earthquake  events  have  been  computed  for  a  network  of 
five  to  ten  linear  arrays.  This  analysis  is  continuing,  testing 
the  correlogram  method  with  known  signals  and  various  noise 
samples  at  selected  signal-noise  ratios. 

E .  Dispersion  Analysis  of  Surface  Waves  from  Deep  Wall  Data 

During  this  reporting  period,  the  digital  programs  for 
computing  deep  well  dispersion  analyses  have  been  improved.  A 
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theoretical  study  of  this  method  is  nearing  completion,  and 
various  analyses  of  deep  well  data  have  been  made  for  other 
investigators . 

F .  Design  of  Optimum  Arrays 

In  the  project,  on  the  optimal  design  of  seismic  arrays, 
the  approach  has  been  to  consider  an  array  as  a  discrete  approxi¬ 
mation  to  an  optimal  wave  number  filter.  This  wave  number  filter 
must  be  designed  from  a  knowledge  of  the  wave  number  power  spectra 
of  signals  and  noise.  Computer  programs  have  been  written  to 
produce  these  spectra.  Studies  are  in  progress  to  determine  the 
best  discrete  approximations.  Recently,  improved  methods  for  ob¬ 
taining  cross  power  spectra  have  been  incorporated  into  the  wave¬ 
number-frequency  spectrum  program. 

G .  Partial  Coherency  Analysis  of  Seismic  Noise 

A  study  is  being  made  of  the  seismic  application  of 
partial  coherencies.  This  study  is  being  carried  out  jointly 
with  the  Measurement  Analsyis  Corporation  under  the  direction  of 
Dr.  Julius  Bendat .  Organizationally,  the  participation  of  MAC 
is  in  the  role  of  consultants.  The  techniques  of  partial  co¬ 
herencies  apply  to  multichannel  filtering  problems.  By  using 
these  methods,  multichannel  transfer  frunctions  may  be  estimated, 
and  tests  are  possible  which  determine  if  the  appropriate  inputs 
and  outputs  have  been  isolated.  The  seismic  applications  will 
include  the  analysis  of  noise  in  arrays,  a  preliminary  study 
has  been  undertaken  which  demonstrates  the  numerical  computations 
involved.  A  report  is  being  prepared  on  the  first  phase  of  the 
preliminary  study.  A  second  report  will  be  concerned  with  the 
actual  numerical  results  obtained  under  controlled  conditions. 
Application  is  simultaneously  being  made  of  these  computational 
tools  to  seismic  data. 

H.  TFSO  Detection  Capability  Study 

Representative  noise  and  seismic  signals  recorded  at 
TFSO  were  combined  as  part  of  a  statistical  study  of  the  VELA 
observatories  detection  capability.  Various  combinations  of 
single  traces  and  summation  traces  are  played  out  on  the  develo- 
corder  and  sent  to  TFO  for  analysis.  A  Geotechnical  Corporation 
program  is  being  used  to  provide  a  statistical  distribution  of 
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the  resulting  measurements  of  ampli  ade  and  period  for  different 
signal  levels. 

W .  SUPPORT  AND  SERVICE  TASKS 

As  part  of  the  contract  work- statement ,  the  SDL  provides 
support,  and  service  functions  for  AFTAC  and  other  VELA  partici¬ 
pants.  With  the  approval  of  AFTAC.  the  SDL  provided  one  or 
more  of  the  following  services  to  organizations  listed  in 
Appendix  A: 

-  copies  of  16  and  35  mm  film 

-  copies  of  existing  composite  analog  tapes 

-  composite  analog  tapes  of  special  events 

-  use  of  1604  computer  for  checking  out  new  programs 
or  running  production  programs 

-  copies  of  digital  programs 

-  digitized  data  in  standard  formats  or  special 
formats  for  use  on  computers  other  than  the  1604 

-  running  SDL  production  programs,  such  as  power 
spectral  density,  and  array  processing  on 
specified  data 

-  digital  x  -  y  plots  of  power  spectra  or  digitized 
data 

Visiting  scientists  have  been  provided  with  space  to  study 
data  and  exchange  information  with  SDL  personnel. 

A.  Data  Library 

The  Data  Library  contains  3400  digitized  seismograms, 

102  digital  computer  programs,  and  192  composite  analog  magnetic 
tapes.  This  data  is  available  for  use  by  the  VELA-UNIFORM  program. 

The  following  additions  were  made  during  this  report 

period: 

1.  Digital  Seismograms  -  1300,  including: 

-  10  Lake  Superior  Explosions 

-  Salmon  long  and  short  period  data 

-  noise  from  several  deep  wells 

-  16  events  recorded  on  British  arrays 


37 


-  Russian  explosions  and  earthquakes 

-  array  data  from  over  150  earthquakes 

2.  Digital  Programs  -  17,  including: 

-  DEPTHMAG  -  This  program  retrieves  epi¬ 
center  information,  from  punched  PDE 
cards,  that  meet  a  certain  criteria 
applied  to  a  number  cf  different  values 
associated  with  each  epicenter.  As  epi¬ 
centers  are  retrieved,  the  program  summa¬ 
rizes  on  magnitude  and  depth. 

-  POINTFLT  ••  This  subroutine  accepts  one 
data  input  at  a  time  and  passes  this  data 
through  a  set  of  recursive  filter  oper¬ 
ations  which  can  be  stated  either  as  a 
time  operation  or  as  a  frequency  operation. 

-  LANDMl  -  Given  the  coefficients  of  a 
numerical  filter  and  consecutive  segments 
of  an  input  data  series,  this  routine 
generates  the  segments  of  the  corresponding 
filtered  series.  The  entire  series  is 
filtered  by  means  of  successive  calls  to 
the  subroutine. 

-  TAPM3RG2  -  To  merge  two  magnetic  tapes, 
made  up  of  epicenter  card  images  on 
magnetic  tap<"  .  The  resulting  merged 
tape  is  in  chronological  sequence  and 
in  BCD  mode  with  end-of-files  between 
months  and  a  double-end-of-f iles  at  the 
end  of  the  year. 

-  COLLATE  -  This  program,  using  a  tape  of 
BCD  card  images  from  two  different  epi¬ 
center  card  lists  that  are  merged  into 
chronological  sequence,  searches  for 
earthquakes  which  are  reported  on  both 
lists.  Matching  of  earthquakes  is  deter¬ 
mined  by  comparison  with  (1)  a  predeter¬ 
mined  origin  time  difference  and  (2)  the 
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sum  of  the  absolute  value  of  the  lati¬ 
tude  difference  and  absolute  value  of  the 
longitude  difference.  Latitude,  Longi¬ 
tude  depth,  magnitude,  and  origin  time 
differences  (residuals)  are  computed 
only  for  matched  pairs  of  earthquakes 
and  are  available  as  output . 

TAPEACT  -  Various  EPICENTER  cards  from 
several  sources  are  written  on  tape  in 
a  standard  format  in  BCD  or  packed  BINARY 
mode .  SEIS  and  GEO  region  numbers  may  be 
added  into  the  format.  Binary  or  BCD 
tapes  may  be  listed  at  end  of  an  update 
run  or  separately.  Packed  binary  tapes 
(formed  by  subroutine  BCDT2)  may  be  con¬ 
verted  back  to  standard  FORTRAN  BCD. 

LAGTIME  -  This  program  combines  a  signal 
and  noise  on  the  same  data  channel  at 
varying  signal-to-noise  ratios.  A  tape 
of  the  new  combination  of  data  is  written 
and  a  plot  tape  is  formed. 

MAKARAY  -  Produces  a  set  of  output  traces 
corresponding  to  any  set  of  element  po¬ 
sitions  desired.  A  noise  background  is 
simulated  by  selecting  ten  traces  from  a 
set  of  twenty  sample  traces,  and  assigning 
to  each  an  arbitrary  velocity,  azimuth, 
and  relative  amplitude.  A  signal  is 
selected  from  a  set,  its  velocity  and 
azimuth  are  chosen  and  a  signal-to-noise 
ratio  is  assigned. 

DISCRETE  -  Expands  either  a  seismic 
signal  or  a  specified  function  in  terms 
of  a  set  of  orthonormalized  exponential 
functions.  Parameters  given  by  the  user 
define  the  orthonormal  functions.  The 
spectrum  of  the  expansion  is  available  as 
an  option. 
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MOD  PV-7  FORTRAN  -  This  program  computes 
for  ail  modes  of  Love  and  Rayleigh  waves 
on  an  elastic  halfspace  of  plane-parallel, 
homogeneous,  isotropic  layers,  the  fallowing: 
phase  velocity,  group  velocity,  and  surface 
orbit  as  functions  of  period  or  frequency; 
amplitude,  the  product  of  vertical  wave 
number  times  layer  thickness,  and  potentials 
as  functions  of  depth;  average  kinetic,  po¬ 
tential,  and  total  energy  densities,  and 
average  horizontal  energy  flux,  corresponding 
to  each  layer;  and  a  summation  of  the  energy 
quantities  from  the  free  surface  .  Provision 
is  made  for  printer  plots  of  phase  and  group 
velocity.  Provision  is  also  made  for  partial 
decoupling  where  amplitude  is  large  m  certain 
channels  at  depth,  relative  to  the  surface 
amplitude . 

LOWPAZ  -  This  program  performs  low  pass  and 
high  pass  filtering  of  data  on  a  SUBSET 
format  tape,  and  produces  another  SUBSET 
tape  and  a  plot  tape  of  the  filtered  data. 

FRAMIS  -  Given  the  parameters  in  the  Laplace 
transform,  to  evaluate  the  amplitude  spectra 
and  phase  of  the  transformation. 

S0RT4  -  To  sort  four  arrays  in  order  of 
increasing  or  decreasing  values  of  one 
of  them. 

XYPLOT  -  To  plot  a  function  of  two  variables 
on  the  printer. 

RESREC  -  Computes  and  plots  the  amplitude 
of  a  recursive  phaseless  digital  filter. 

FLAP  -As1'  jrdinate  control  system  which 
allows  chaining  of  programs  in  any  com¬ 
bination  to  perform  a  specified  analysis. 
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-  PHASDIST  -  Generates  phasing  information 
for  an  arbitrary  array  of  je ismometers , 
which  may  be  input  to  the  computer 
program,  ARRAY. 


Modifications  were  made  to  the  ARRAY  program  in  order 
to  process  array  data  using  the  British  filtering  and  smoothing 
techniques.  In  addition,  modifications  were  made  to  the  following 
programs : 

1.  The  analysis  of  variance  program  to  compute  and 
print  . 

2.  Prediction  error  fill,  ring  program  to  output  results 
in  a  form  which  can  then  be  processed  with  the 
analysis  of  variance  program. 

3.  Program  DEPTHMAG  to  work  with  different  types  •„  * 
input  data,  an  to  search  -nd  r-' '  neve  on  the 
parameter  which  specifies  the  source  of  the 
magnitude  calculation. 

3.  Analog  Composite  Tapes  -  23,  including: 

a.  Made  by  Seismic  Data  Laboratory 


-  SALMON 

-  SS  VILLAGE 

-  Russian  explosion 

-  Special  composites  made  for  the  U.  S. 
Coast  &  Geodetic  Survey,  Roland  F.  Beers, 
Inc.,  Texas  Instruments,  Inc.,  and  i'.amont 
Geophysical  Observatory 

b.  Made  by  the  Geotechnical  Corporation 


-  YUBA 

-  SATSOP 

-  KINNEBEC 

-  GERBIL 

-  HUTTA 

-  STONES 

-  FERRET 

-  MATACO 


-  KOHOCTON 

-  PECAN 

-  CLEARWATER 

-  ANCHOVY 

-  SARDINE 

-  EAGLE 

-  GRUNION 

-  FORE 
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4 .  Other  D ata  Available  to  VELA  Participants 


The  following  data  is  stored  on  magnetic  tape: 

-  U.  S.  Coast  and  Geodetic  Survey  epi- 
enters  from  I960  to  date 

-  Earthquake  Bulletin  Data  from  the 
LRSM  stations  (January  1962  to  October 
1964) 

-  Earthquake  Bulletin  Data  from  the  VELA 
observatories  (February  1963  to  Novembe 
1964) 

-  Amplitude,  period,  and  arrival  times  of 
phases  from  50  U.  S.  shots 

B .  Digital  Computer  Operations 

The  use  of  digital  computer  machine  data  processing 
techniques  for  in-houre  projects  and  for  other  VELA  participants 
has  been  increasing,  ana  recently  a  second  shift  digital  compute 
operation  was  started  to  handle  this  increased  work-load. 

The  SDL  took  over  the  maintenance  of  the  1604  digital 
computer  system  in  January  1965,  after  a  six  month  period  of  on- 
t he- job  training.  Necessary  digital  spare  parts  were  purchased 
and  a  punched  card  inventory  system  has  been  set  up  for  all 
computer  spare  parts. 

C .  Data  Compression 

This  is  continuing  as  a  routine  operation,  and 
production  is  maintained  at  the  level  needed  to  meet  the 
requirements  of  the  field  operations  (LRSM  teams  and  U.  S. 
observatories)  and  the  Seismic  Data  Laboratory. 

D .  Equipment  Modifications 

Six  additional  Dynamic  amplifiers  were  installed  in 
the  playback  system  to  give  SDL  the  capability  of  making  12 
channel  playouts  with  the  CEC  oscillograph. 
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The  air  conditioning  system  in  the  computer  room 
was  modified  so  as  to  be  synchronized  with  the  operation  of  the 
digital  computer. 

A  cycle  adder  i being  installed  on  the  Sangamo  tape 
transports.  This  device  v.ill  permit  the  SDL  to  synchronize 
recording  times  when  processing  data  simultaneously  from  two 
tape  recorders  . 
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ORGANIZATIONS  RECEIVING  SDL  DATA  SERVICES 


OCTOBER  1964  -  MARCH  1965 


Albuquerque  Seismological  Center 
Bolt,  Beranek  &  Newman,  Inc. 
Bureau  of  Standards 
California  Institue  of  Tech. 
Canadian  Dept,  of  Mines 
Engineering  Physics 
General  Atronics  Corp. 
Geotechnical  Corp. 

Institue  for  Meteorology  and 
Geophysics  (Germany) 
International  Seismological 
Research  Centre  (Scotland) 
Lamont  Geological  Observatory 
Pennsylvania  State  University 

Oregon  State 


Princeton  University 
Roland  F.  Beers  Co. 
Southwest  Center  for 
Advanced  Studies 
St.  Louis  University 
Texas  Instrument,  Inc. 
Underwater  Systems,  Inc. 
University  of  California 
University  of  Michigan 
University  of  Sweden 
U.  S.  Coast  &  Geodetic 
Survey 

Vitro  Corporation 
Xavier  University 
Univeristy 
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detection;  25  February  1965. 
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Seismic  waves  from  the  SS  VILLAGE  explosion; 
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Texas-Louisiana  earthquake;  9  Felruary  196a. 
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Van  Nostrand, 


R.;  Reverberation  effects  on  signal 
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Booker,  A. 
Booker,  A. 
Booker,  A. 

Claerbout , 

Dean,  W.  C 


H.;  A  statistical  discriminator?  22  November  1962. 

H.;  Estimation  of  network  capability;  20  January  1964. 

H.;  Parameter  dependence  of  integrated  seismic  phases; 
25  May  1964. 

J.  F.?  Detection  of  P  waves  from  weak  sources  at 
great  distances;  13  September  1963. 

;  Inverse  filtering  of  seismic  signals?  1  August 
1962  , 


Dean,  W .  C.;  Cor.velograms  of  Mississippi  from  orthogonal 

filters;  1  March  1963. 

FI inn,  E.  A.;  Confidence  regions  and  error  determinations  for 

seismic  event  lo <  tion;  24  April  1964. 

Flinn,  E.  A.  &.  Engdahl,  E.  R.*?  A  proposed  basis  for  geographical 

and  seismic  regionalization?  9  March  1964. 


Martinek,  J.,  and  Thielman,  H.  P.;  Special  sphere  and  circle 

theorems  for  the  general  biharmonic  field; 

8  July  1963. 


Martinek,  J.  &  Thielman,  H.  P.;  Special  sphere  and  circle 

theorems  for  the  general  biharmonic  equation  - 
interior  and  exterior  problems;  11  December  1963. 


Martinek,  J.  &  Thielman,  H.  P.;  Propagation  of  waves  from  a 

spherical  cavity  into  concentric  spheres? 

13  March  1964. 

Martinek,  J.,  and  Thielman,  H.  P.;  Methods  for  obtaining  new 

solutions  to  the  harmonic  and  wave  equations; 
30  April  1964. 

Mims,  C.  H.?  A  sum-of-squares  method  of  seismic  phase 

identification;  27  November  1962. 
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Mims <  C.  H.;  Detectability  of  first  motion;  15  August  1963. 

Mitronovas,  W.;  Identification  of  nuclear  explosions  using 

ratios  of  smoothed  and  combined  seismic  phases; 

28  June  1963 . 

Mitronovas,  W.;  Identification  of  nuclear  explosions  on  the 

shape  of  initial  P-phase  for  teleseismic  events; 

28  August  1964. 

Pilotte,  F.  F.*;  and  Flinn,  E.  A.;  A  preliminary  evaluation 

of  a  method  for  linear  array  processing; 

31  January  1964. 

Sax,  R.  L.;  Attenuation  of  head  waves  as  a  function  of 

frequency  and  distance;  28  February  1963. 

Sax,  R.  L.,  and  Hartenberger ,  R.  A.;  Theoretical  estimation  of 

the  signal -to-noise  ratio  as  a  function  of  depth; 

3  April  1963. 

Sax,  R.  L.  and  Hartenberger,  R.  A.:  Theoretical  prediction  of 

seismic  noise  in  a  deep  borehole;  24  September  1963. 

Sax,  R.  L.,  and  Hartenberger,  R.  A.;  Seismic  noise  attenuation 

in  unconsolidated  material;  28  August  1964. 

Thielman,  H.  P.,  and  Sax,  R.  L.;  Relation  giving  the  phase 

velocity  in  terms  of  the  group  velocity; 

28  September  1962. 

Van  Nostrand,  R.;  Synthetic  earthquake  seismograms  for  tele- 

seismic  distances;  14  May  1964. 

Van  Nostrand,  and  Helterbrand,  W.*;  A  comparative  study  of 

the  SHOAL  event;  29  August  1964. 
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